Modelos indivuais estocásticos estruturados

A importância da estrutura da rede de conexões

• Lembremos da figura a seguir, retirada do livro Kiss, Miller & Simon (2017), que mostra a importância da estrutura da rede conexões na evolução da epidemia.

• O modelo dele não é individual, mas é um compartimental estruturado estocástico, com muitos compartimentos ($100$ mil nós), que carrega bastante estrutura.

Image(filename=path.join('..', 'input', 'networks', 'SIR_on_different_networks.jpg'), width=800)

Exemplos de simulações individuais

• A seguir, três vídeos retirado das informações suplementares do artigo de Ferguson et al. (2005), sobre estratégias para conter pandemias de influenza na sudeste asiático.

• Observem que eles já estavam pensando nesses estratégias!

• Vieram motivadas pela epidemia de SARS que apareceu no final de 2002 e só foi considerada contida em meados de 2003, com vários casos ocorrendo até o iníco de 2004.

• Simulação com 85 milhões de indivíduos!

• A primeira, sem estratégia de contenção da epidemia.

• A segunda, estratégia mais extensa, bem sucedida.

• A terceira, estratégia em apenas um país, mal sucedida.

Sem estratégia

• $R_0=1.5$.

• Video mostra 100 dias de evolução.

• Cores:

  • Cinza: densidade de suscetíveis em escala logarítmica
  • Vermelho: áreas infectadas
  • Verde: áreas recuperadas
  • Azul: áreas em tratamento

show_video(path.join('..', 'input', 'networks', '41586_2005_BFnature04017_MOESM2_ESM.mp4'), '50%')

Estratégia bem sucedida

• $R_0=1.8$;

• Estratégia implementada na Tailândia e em parte dos países vizinhaça, dentro de um raio de 100km da Tailândia.

• Distanciamento social

• Quarentena de casos detectados

• Profilaxia em um raio de 5km radial de casos dectados

• Fechamento de 90% de escolas e de 50% dos lugares de trabalho, em um raio de 5km.

• Video mostra 300 dias de evolução.

• Cores:

  • Cinza: densidade de suscetíveis em escala logarítmica
  • Vermelho: áreas infectadas
  • Verde: áreas recuperadas
  • Azul: áreas em tratamento

show_video(path.join('..', 'input', 'networks', '41586_2005_BFnature04017_MOESM3_ESM.mp4'), '80%')

Estratégia mal-sucedida

• $R_0=1.7$.

• Estratégias implementadas em apenas um país: Tailândia

• Distanciamento social

• Quarentena de casos detectados

• Profilaxia em um raio de 5km radial de casos dectados

• Fechamento de 90% de escolas e 50% de lugares de trabalho no raio de 5km.

• Video mostra 300 dias de evolução.

• Cores:

  • Cinza: densidade de suscetíveis em escala logarítmica
  • Vermelho: áreas infectadas
  • Verde: áreas recuperadas
  • Azul: áreas em tratamento

show_video(path.join('..', 'input', 'networks', '41586_2005_BFnature04017_MOESM4_ESM.mp4'), '50%')

Modelagem

Parâmetros

Nas modelagens que faremos a seguir:

• População: $300$ indivíduos;

• $\beta= 0.5$ (taxa diária de contágio);

• $\gamma=0.1$ (taxa de recuperação - 10 dias);

• $R_0 = 5$ (número básico de reprodução).

• $I_0 = 8$ (número inicial de infectados);

• $T = 80$ (80 dias de simulação);

Nos exemplos com rede, usamos

• $\kappa = $ número médio de conexões por indivíduo (arestas por vértice);

• $\mu = \beta/\kappa$ (contágio por contato).

num_pop = 300
beta = 0.5
gamma = 0.1
R0 = beta/gamma
I_0 = 8
T = 80
tempos = np.array(range(T+1))
num_sim = 60
print(f'Número básico de reprodução = {R0:.2f}')

Número básico de reprodução = 5.00

Estado da população

• Temos um índice $i=1, \ldots, N$ para cada indivíduo.

• O estado do sistema todo é representado por um vetor $\mathbf{x} = (x_1, \ldots, x_N)$ de dimensão $N$.

• Cada coordenada $x_i$ do vetor $\mathbf{x}$ representa o estado do individuo $i$.

• Cada indivíduo pode estar em um dos possívei estágios: suscetível, infectado ou removido..

• Representamos isso permitindo que $x_i$ assuma três valores, cada um correspondendo a um dos estágios:

  • $x_i=1:$ indivíduo suscetível;

  • $x_i=2:$ indivíduo infectado;

  • $x_i=3:$ indivíduo removido. O estado da população é um vetor

# estado inicial da população
pop_0 = np.ones(num_pop)
infectados_0 = np.random.choice(num_pop, I_0)
pop_0[infectados_0] = 2*np.ones(I_0)
attr_pop_0 = dict([(i, {'estado': int(pop_0[i])}) for i in range(num_pop)])
print(pop_0)

[1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 2. 1. 1. 1. 1. 1. 1. 1. 1. 1.
 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 2. 1. 2. 1. 1. 1. 1. 1. 1.
 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 2. 1. 1. 1.
 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
 1. 1. 1. 2. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 2. 1. 1. 1. 1. 1. 1.
 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 2. 1. 1. 1. 1. 1. 1.
 1. 1. 1. 1. 1. 1. 1. 1. 1. 2. 1. 1.]

Evolução temporal

• Método $\tau$-leap ($\tau$-leap method):

  • Evolução discreta no tempo.

  • Múltiplos eventos ocorrendo a cada passo de tempo.

  • Taxas de ocorrência de algum evento se tornam probabilidades da forma

    $$ \mathbb{P}(\textrm{evento}) = 1 - e^{\displaystyle - \textrm{taxa} \,\times \,\textrm{passo de tempo}}. $$

  • Isto segue de $y' = -\alpha y$ $\Longrightarrow$ $y(t+\delta t) = y(t) e^{-\alpha \delta t}$ $\Longrightarrow$ $\Delta y = y(t+\delta t) - y(t) = (1 - e^{-\alpha \delta t})y(t)$.

• Passos muito grandes diminuem a acurácia do método

• Passos muito curtos o tornam bastante custoso caso haja um número muito grande de indivíduos

• Veja Box 7.3, página 270, de Keeling & Rohani (2007).

• As simulações acima na Tailândia foram com passos de 1/4 de dia, ou 6 horas.

• Vamos usar passos de um dia apenas para efeito de ilustração.

Funções de evolução e de análise

Definimos algumas funções para facilitar as simulações.

• evolucao_SIR

• analise_grafo

• evolucao_grafo

def evolucao_SIR(pop_0, beta, gamma, tempos, plot=False):

    SIR_Compartimental = namedtuple('SIR_Compartimental', 
                                 [
                                     'pop_0',
                                     'num_pop',
                                     'beta',
                                     'gamma',
                                     't_instantes',
                                     'S',
                                     'I', 
                                     'R',
                                     'X'
                                 ])        
          
    def diferencial(t, X, N, beta, gamma):
        S, I = X
        dXdt = [- beta*I*S/N, beta*I*S/N - gamma*I]
        return dXdt        
    
    num_pop = sum(pop_0)
    y0 = pop_0[0:2]
    sol = solve_ivp(diferencial, t_span=[tempos[0],tempos[-1]], 
                    y0 = y0, t_eval = tempos,
                    args=(num_pop, beta, gamma))
    
    resultado = SIR_Compartimental(pop_0, num_pop, beta, gamma, tempos,
                                   sol.y[0], sol.y[1], num_pop - sol.y[0] - sol.y[1], sol)
    return resultado

def analise_grafo(G, info=True, node_size=0, pos=None, hist=False):

    num_vertices = G.number_of_nodes()
    num_arestas = G.number_of_edges()
    num_arestas_por_vertice = [j for i,j in G.degree]
    num_medio_conexoes = 2*num_arestas/num_vertices

    if info:
        print(f'Número de arestas: {num_arestas}')
        print(f'Número de vértices: {num_vertices}')
        print(f'Número médio de conexões por indivíduo: {num_medio_conexoes:.1f}')

    if node_size:
        color = ['tab:blue', 'tab:red', 'tab:green']
        pop_estado = nx.get_node_attributes(G,'estado')
        color_map = [color[pop_estado[j]-1] for j in range(num_vertices)]
        plt.figure(figsize=(10,6))
        if pos:
            nx.draw(G, pos, node_size=node_size, node_color=color_map, alpha=0.5)
        else:
            nx.draw(G, node_size=node_size, node_color=color_map, alpha=0.5)
        plt.title('Rede de indivíduos e de suas conexões', fontsize=16)
        plt.show()

    if hist:
        plt.figure(figsize=(10,6))
        plt.hist(num_arestas_por_vertice, 50, facecolor='g', alpha=0.75)
        plt.xlabel('num. arestas', fontsize=14)
        plt.ylabel('num. vertices', fontsize=14)
        plt.title('Histograma do número vertices por número de arestas', fontsize=16)
        plt.show()
    return num_medio_conexoes

def evolucao_grafo_estruturado(pop_0, mu, gamma, G, tempos, num_sim, show=''):
    """Evolução temporal da epidemia em um grafo estruturado.


    Entrada:
        pop_0: numpy.array
            state of the population, with 
                1: suscetível
                2: infectado
                3: recuperado ou removido

        mu: float
            probabilidade de contágio por contato

        gamma: float
            taxa de recuperação por unidade de tempo

        G: numpy.array
            grafo de conexões

        tempos: numpy.array
            instantes de tempo

        num_sim: int
            número de simulações

        show: str
            indica se é para exibir um gráfico e de que tipo:
                - 'nuvem': exibe uma nuvem com todas as simulações e o valor médio em destaque
                - 'sd': exibe o valor médio com um intervalo de confiança dado pelo desvio padrão
                - 'medio': exibe apenas o valor médio
                - '': não exibe gráfico algum.

    Saída
        X: class.SIR_Individual
            Uma instância da classe `SIR_Individual` com os seguintes atributos:
                pop_0:
                num_pop:
                mu:
                gamma:
                G:
                tempos:
                num_sim:
                S_mean:
                I_mean: 
                R_mean:
                I_sigma:
                R_sigma:
                S_sigma:
    """

    def passo_grafo_estruturado(num_pop, populacao, A_adj, prob_infeccao, prob_recuperacao):
        """
        DEFINITIVAMENTE NÃO ESTÁ OTIMIZADO PARA VELOCIDADE DE PROCESSAMENTO E MEMÓRIA
        
        Foi otimizado para facilidade de construção.
        
        Por exemplo, geramos muito mais números aleatórios do que realmente precisamos. 
        Precisamos de, no máximo, num_edges + num_pop número aleatórios, mas geramos num_pop**2.
        """
        
        # gera uma matriz cheia aleatória (números em [0.0,1.0))
        A_random = np.random.rand(num_pop, num_pop)
        
        # filtra de acordo com as conexões
        A_graph_random = np.multiply(A_adj, A_random)
        
        # copia a população (acho que é desnecessário)
        pop_aux = np.copy(populacao)
        
        # separa os suscetíveis, criando um vetor de 1's e 0's, se for, ou não, suscetível
        pop_suscetiveis = np.select([pop_aux==1], [pop_aux])
        
        # separa os infectados, criando um vetor de 1's e 0's, se for, ou não, infectável
        pop_infectados = np.select([pop_aux==2], [pop_aux])/2
        
        # cria uma matriz de risco, mantendo apenas as conexões que envolvem um infectado
        A_risco = np.multiply(np.tile(pop_infectados, (num_pop, 1)), A_graph_random)
        
        # cria matriz de contatos entre um suscetível e um infectado
        A_contatos = np.multiply(np.tile(pop_suscetiveis, (num_pop, 1)).transpose(), A_risco)
        
        # cria uma matriz de 1's e 0's, indicando se houve, ou não, contágio
        A_infectados = np.select([A_contatos > 1-prob_infeccao], [np.ones([num_pop, num_pop])])

        # obtém novos infectados
        pop_novos_infectados = np.select([np.sum(A_infectados, axis=1)>0], [np.ones(num_pop)])
        
        # filtra matriz aleatória com a diagonal
        pop_recuperando = pop_infectados @ np.multiply(np.eye(num_pop), A_random)
        
        # obtém novos recuperados
        pop_novos_recuperados = np.select([pop_recuperando > 1-prob_recuperacao], [np.ones(num_pop)])
        
        # atualiza população adicionando um aos que avançaram de estágio
        populacao_nova = pop_aux + pop_novos_infectados + pop_novos_recuperados

        # Observe que cada elemento da matriz aleatória é usado apenas uma vez, garantindo
        # a independência desses eventos aleatórios (tanto quanto se leve em consideração
        # que os números gerados são pseudo-aleatórios)

        return populacao_nova

    
    # confere se escolha para `show` é válida
    if show:
        assert(show in ('nuvem', 'sd', 'media')), 'Valor inválido para argumento `show`.'

    # atributos de saída
    SIR_Individual = namedtuple('SIR_Individual', 
                                [
                                    'pop_0',
                                    'num_pop',
                                    'mu',
                                    'gamma',
                                    'G',
                                    'tempos',
                                    'num_sim',
                                    'S_mean',
                                    'I_mean', 
                                    'R_mean',
                                    'S_sigma',
                                    'I_sigma', 
                                    'R_sigma'
                                ])
    
    # número de indivíduos da população
    num_pop = len(pop_0)
    I_0 = np.count_nonzero(pop_0==2)        

    # número de instantes no tempo e passos de tempo 
    num_t = len(tempos)
    passos_de_tempo = tempos[1:] - tempos[:-1]

    # inicializa variáveis para o cálculo da média
    S_mean = np.zeros(num_t)
    I_mean = np.zeros(num_t)
    R_mean = np.zeros(num_t)

    # inicializa variáveis para o cálculo do desvio padrão
    S_sigma = np.zeros(num_t)
    I_sigma = np.zeros(num_t)
    R_sigma = np.zeros(num_t)
    
    # obtém matriz de adjacências e o número médio de vizinhos do grafo
    A_adj = nx.to_numpy_matrix(G)

    # prepara gráfico se necessário
    if show in ('nuvem', 'sd', 'media'):    
        # inicializa figura e define eixo vertical
        plt.figure(figsize=(12,6))
        plt.ylim(0, num_pop)
        plt.xlim(tempos[0], tempos[-1])
    
    if show == 'nuvem':
        # alpha para a nuvem de gráficos das diversas simulaçõe
        alpha = min(0.2, 5/num_sim)    

    # simulações
    for j in range(num_sim):

        # inicializa população de cada simulação
        populacao = np.copy(pop_0)
        S = np.array([num_pop - I_0])
        I = np.array([I_0])
        R = np.array([0])
     
        
        # evolui o dia e armazena as novas contagens
        for dt in passos_de_tempo:

            # calcula probabilidades
            prob_infeccao = 1 - np.exp(-mu*dt)
            prob_recuperacao = 1 - np.exp(-gamma*dt)

            populacao = passo_grafo_estruturado(num_pop, populacao, A_adj,
                                                prob_infeccao, prob_recuperacao)
            S = np.hstack([S, np.count_nonzero(populacao==1)])
            I = np.hstack([I, np.count_nonzero(populacao==2)])
            R = np.hstack([R, np.count_nonzero(populacao==3)])
            
        # adiciona as contagens dessa simulação para o cálculo final da média
        S_mean += S
        I_mean += I
        R_mean += R

        # adiciona as contagens dessa simulação para o cálculo final do desvio padrão
        S_sigma += S ** 2
        I_sigma += I ** 2
        R_sigma += R ** 2

        if show == 'nuvem':
            # exibe os gráficos dos dados de cada simulação
            plt.plot(tempos, S, '-', color='tab:green', alpha=alpha)
            plt.plot(tempos, I, color='tab:red', alpha=alpha)
            plt.plot(tempos, R, '-', color='tab:blue', alpha=alpha)
            plt.plot(tempos, num_pop - S, '-', color='tab:gray', alpha=alpha)

    # divide pelo número de evoluções para obter a média
    S_mean /= num_sim
    I_mean /= num_sim
    R_mean /= num_sim

    # ajusta o calcula do desvio padrão
    S_sigma = ( S_sigma / num_sim - S_mean**2 )**.5
    I_sigma = ( I_sigma / num_sim - I_mean**2 )**.5
    R_sigma = ( R_sigma / num_sim - R_mean**2 )**.5

    # exibe os gráficos das médias
    if show in ('nuvem', 'sd', 'media'):
        plt.plot(tempos, S_mean, '-', color='tab:green', label='suscetíveis')
        plt.plot(tempos, I_mean, '-', color='tab:red', label='infectados')
        plt.plot(tempos, R_mean, '-', color='tab:blue', label='recuperados')
        plt.plot(tempos, num_pop - S_mean, '-', color='tab:gray', label='inf.+ rec.')
        
    if show == 'sd':
        plt.fill_between(tempos, S_mean - S_sigma, S_mean + S_sigma, facecolor='tab:green', alpha = 0.2)
        plt.fill_between(tempos, I_mean - I_sigma, I_mean + I_sigma, facecolor='tab:red', alpha = 0.2)
        plt.fill_between(tempos, R_mean - R_sigma, R_mean + R_sigma, facecolor='tab:blue', alpha = 0.2)
        plt.fill_between(tempos, num_pop - S_mean - S_sigma, num_pop - S_mean + S_sigma, facecolor='tab:gray', alpha = 0.2)

        
    if show == 'nuvem':
        plt.title('Evolução do conjunto de simulações e da média', fontsize=16)
    elif show == 'sd':
        plt.title('Evolução da média, com o desvio padrão', fontsize=16)
    elif show == 'media':
        plt.title('Evolução da média das simulações', fontsize=16)

        
    # informações para o gráfico
    if show in ('nuvem', 'sd', 'media'):
        plt.xlabel('tempo', fontsize=14)
        plt.ylabel('número de indivíduos', fontsize=14)
        plt.legend(fontsize=12)
        plt.show() 

    resultado = SIR_Individual(
        pop_0,
        num_pop,
        mu,
        gamma,
        G,
        tempos,
        num_sim,
        S_mean,
        I_mean, 
        R_mean,
        I_sigma,
        R_sigma,
        S_sigma
    )

    return resultado

SIR Compartimental

Começamos com a simulação SIR Compartimental, para servir de comparação com as simulações indi

X = evolucao_SIR(pop_0 = [num_pop-I_0, I_0, 0], beta=beta, gamma=gamma, tempos = tempos)

plt.figure(figsize=(12,6))

# exibe os gráficos das médias
plt.plot(tempos, X.S, 'tab:blue', label='SIR: suscetíveis')
plt.plot(tempos, X.I, 'tab:red', label='SIR: infectados')
plt.plot(tempos, X.R, 'tab:green', label='SIR: recuperados')
plt.plot(tempos, X.I + X.R, 'tab:gray', label='SIR: inf. + sus.')

# informações para o gráfico
plt.xlabel('dias', fontsize=14)
plt.ylabel('número de indivíduos', fontsize=14)
plt.title('Evolução da população em cada compartimento da simulação do modelo SIR', fontsize=16)
plt.legend(loc='best', fontsize=12)
plt.show() 

Redes

• Agora vamos fazer simulações em algumas redes.

• Neste primeiro momento, consideramos apenas redes fixas e com contatos uniformes.

• Posteriormente, devemos incluir outras conexões e pesos

  • redes para as conexões frequentes
    • rede estruturada fixa com arestas com pesos determinados pela relação entre os indivíduos
    • Exemplos:
      • familiares, com um peso maior;
      • trabalho e escola;
      • amigos;
  • redes para as conexões aleatórias
    • Exemplos:
      • praças, shoppings e parques
      • em meios de transporte;
    • rede completa geométrica, com pesos dependendo da distância;

• Para a criação das redes, usaremos funções descritas em NetworkX 2.4: Graph Generators para gerar alguns grafos sintéticos.

Simulação em uma rede regular

• Começamos com uma rede regular.

• Cada indivíduo (vértice) tem o mesmo número de conexões (arestas).

• Usamos a função nx.gnp_random_graph.

Criando a rede e analisando seus dados

G_r = nx.random_regular_graph(d=num_pop-1, n=num_pop, seed=12)
nx.set_node_attributes(G_r, attr_pop_0)

num_medio_conexoes_r = analise_grafo(G_r, info=True, node_size=80, pos=None, hist=True)

Número de arestas: 44850
Número de vértices: 300
Número médio de conexões por indivíduo: 299.0

Simulação

%time
X_r = evolucao_grafo_estruturado(pop_0, beta/num_medio_conexoes_r, gamma, G_r, tempos,
                                 num_sim, show='nuvem')

CPU times: user 5 µs, sys: 1 µs, total: 6 µs
Wall time: 11 µs

Comparação com o SIR compartimental

plt.figure(figsize=(12,6))
plt.ylim(0, num_pop)

# exibe os gráficos das médias
plt.plot(tempos, X_r.S_mean, '.', color='tab:green', label='média suscetíveis')
plt.plot(tempos, X_r.I_mean, '.', color='tab:red', label='média infectados')
plt.plot(tempos, X_r.R_mean, '.', color='tab:blue', label='média recuperados')

plt.plot(tempos, X.S, 'tab:green', label='compart. suscetíveis')
plt.plot(tempos, X.I, 'tab:red', label='compart. infectados')
plt.plot(tempos, X.R, 'tab:blue', label='compart. recuperados')

# informações para o gráfico
plt.xlabel('dias', fontsize=14)
plt.ylabel('número de indivíduos', fontsize=14)
plt.title(f'Evolução do SIR compartimental e da média do SIR individual estocástico em rede regular', fontsize=16)
plt.legend(loc='best', fontsize=12)
plt.show() 

Simulação em uma rede do tipo Ërdos-Renyi

• Rede aleatória com a propriedade do número de conexões de cada indivíduo (arestas para cada vértice) ser próximo do número médio.

• A função erdos_renyi_graph() do módulo networkx cria aleatoriamente uma rede desse tipo.

• Recebe como argumentos o número de vértices e uma dada densidade média do número de arestas.

• A função to_numpy_matrix do networkx retorna a matriz de adjacências, com as conexões (arestas) entre os indivíduos (vértices).

Criando a rede e analisando seus dados

G_er = nx.erdos_renyi_graph(num_pop, 6/num_pop, seed=12)
nx.set_node_attributes(G_er, attr_pop_0)

num_medio_conexoes_er = analise_grafo(G_er, info=True, node_size=80, pos=None, hist=True)

Número de arestas: 821
Número de vértices: 300
Número médio de conexões por indivíduo: 5.5

Simulação

%time 
X_er = evolucao_grafo_estruturado(pop_0, beta/num_medio_conexoes_er, gamma, G_er, tempos,
                                  num_sim, show='sd')

CPU times: user 3 µs, sys: 1e+03 ns, total: 4 µs
Wall time: 10 µs

Comparação com o SIR compartimental

plt.figure(figsize=(12,6))
plt.ylim(0, num_pop)

# exibe os gráficos das médias
plt.plot(tempos, X_er.S_mean, '.', color='tab:green', label='média suscetíveis')
plt.plot(tempos, X_er.I_mean, '.', color='tab:red', label='média infectados')
plt.plot(tempos, X_er.R_mean, '.', color='tab:blue', label='média recuperados')

plt.plot(tempos, X.S, 'tab:green', label='compart. suscetíveis')
plt.plot(tempos, X.I, 'tab:red', label='compart. infectados')
plt.plot(tempos, X.R, 'tab:blue', label='compart. recuperados')

# informações para o gráfico
plt.xlabel('dias', fontsize=14)
plt.ylabel('número de indivíduos', fontsize=14)
plt.title('Evolução do SIR compartimental\n e da média do SIR individual estocástico em rede Ërdos-Renyi', fontsize=16)
plt.legend(loc='best', fontsize=12)
plt.show() 

Simulação com um grafo Barabasi-Albert

• Possui indivíduos com muitas conexões.

G_ba = nx.barabasi_albert_graph(num_pop, 6, seed=1)
nx.set_node_attributes(G_ba, attr_pop_0)
num_medio_conexoes_ba = analise_grafo(G_ba, info=True, node_size=100, pos=None, hist=True)

Número de arestas: 1764
Número de vértices: 300
Número médio de conexões por indivíduo: 11.8

Simulação

X_ba = evolucao_grafo_estruturado(pop_0, beta/num_medio_conexoes_ba, gamma,G_ba, tempos,
                                  num_sim, show='sd')

Comparação com o SIR Compartimental

plt.figure(figsize=(12,6))

# exibe os gráficos das médias
plt.plot(tempos, X_r.I_mean, '.', label='Rede regular: média infectados')
plt.plot(tempos, X_er.I_mean, '.', label='Ërdos-Remyi: média infectados')
plt.plot(tempos, X_ba.I_mean, ':', label='Barabasi-Albert: média infectados')
plt.plot(tempos, X.I, 'tab:red', label='SIR: infectados')

# informações para o gráfico
plt.xlabel('dias', fontsize=16)
plt.ylabel('número de indivíduos', fontsize=16)
plt.title('Evolução média do SIR compartimental e do SIR individual estocástico\n'
          + f'População {num_pop}; Infectados inicialmene: {I_0}\n'
          + f'beta = {beta:.3f}; gamma = {gamma:.3f}; $R_0$ = {beta/gamma:.3f}', fontsize=18)
plt.legend(loc='best', fontsize=12)
plt.show() 

Simulação com um grafo Barbell

• Rede com duas subredes regulares completas e conectadas por uma cadeia simples.

G_b = nx.barbell_graph(int(num_pop/2)-2, 4)
nx.set_node_attributes(G_b, attr_pop_0)
num_medio_conexoes_b = analise_grafo(G_b, info=True, node_size=30, pos=None, hist=True)

Número de arestas: 21761
Número de vértices: 300
Número médio de conexões por indivíduo: 145.1

Simulação

X_b = evolucao_grafo_estruturado(pop_0, beta/num_medio_conexoes_b, gamma, G_b, tempos,
                                 num_sim, show='sd')

Comparação com o SIR Compartimental e grafos anteriores

plt.figure(figsize=(12,6))
#plt.ylim(0, num_pop/2)

# exibe os gráficos das médias
plt.plot(tempos, X_er.I_mean, '.', label='Ërdos-Remyi: média infectados')
plt.plot(tempos, X_ba.I_mean, ':', label='Barabasi-Albert: média infectados')
plt.plot(tempos, X_b.I_mean, '.-', label='Barbell: média infectados')
plt.plot(tempos, X.I, 'tab:red', label='SIR: infectados')

# informações para o gráfico
plt.xlabel('dias', fontsize=14)
plt.ylabel('número de indivíduos', fontsize=14)
plt.title('Evolução do SIR compartimental e da média do SIR individual estocástico em diferentes redes\n'
          + f'População {num_pop}; Infectados inicialmene: {I_0}\n'
          + f'beta = {beta:.3f}; gamma = {gamma:.3f}; $R_0$ = {beta/gamma:.3f}', fontsize=16)
plt.legend(loc='best', fontsize=12)
plt.show() 

Simulação com um grafo geométrico

• Vamos considerar, agora, um grafo geométrico aleatório construído via função random_geometric_graph (veja, também, o seguinte exemplo).

G_g = nx.random_geometric_graph(num_pop, 0.105, seed=2024)
nx.set_node_attributes(G_g, attr_pop_0)
pos_g = nx.get_node_attributes(G_g, "pos")
num_medio_conexoes_g = analise_grafo(G_g, info=True, node_size=80, pos=pos_g, hist=True)

Número de arestas: 1401
Número de vértices: 300
Número médio de conexões por indivíduo: 9.3

Simulação

X_g = evolucao_grafo_estruturado(pop_0, beta/num_medio_conexoes_g, gamma, G_g, tempos,
                                 num_sim, show='sd')

Comparação com o SIR Compartimental e os grafos anteriores

plt.figure(figsize=(12,6))
#plt.ylim(0, num_pop/2)

# exibe os gráficos das médias
plt.plot(tempos, X_er.I_mean, '.', label='Ërdos-Remyi: média infectados')
plt.plot(tempos, X_ba.I_mean, ':', label='Barabasi-Albert: média infectados')
plt.plot(tempos, X_b.I_mean, '.-', label='Barbell: média infectados')
plt.plot(tempos, X_g.I_mean, 'o-', label='Geométrico: média infectados')
plt.plot(tempos, X.I, 'red', label='SIR: I')

# informações para o gráfico
plt.xlabel('dias', fontsize=14)
plt.ylabel('número de indivíduos', fontsize=14)
plt.title('Evolução do SIR compartimental clássico e da média do SIR individual estocástico em diferentes redes\n'
          + f'População {num_pop}; Infectados inicialmene: {I_0}\n'
          + f'beta = {beta:.3f}; gamma = {gamma:.3f}; $R_0$ = {beta/gamma:.3f}', fontsize=16)
plt.legend(loc='best', fontsize=12)
plt.show() 

A fazer

• Separar as conexões fixas de acordo com o tipo (núcleo familiar, trabalho/escola, local, etc.)

• Incluir pesos para as conexões fixas.

• Incluir um grafo completo de conexões com peso dado por um núcleo em função da distância e com decaimento em lei de potência.

• Simular o efeito de fechamento de estabelecimentos, distânciamento social e métodos de prevenção nos parâmetros ligados a cada tipo de conexão.

• Obter um grafo mais realista para a cidade do Rio de Janeiro.

Trabalho 7

Questões:

1. Adaptar o modelo para incluir os estados SEIR.

2. Comparar esse modelo com o SEIR compartimental.

3. Adaptar o modelo para incluir os estados SEIAR.

4. Comparar esse modelo com o SEIAR compartimental.

5. Adaptar o modelo para incluir os estados SEIR-QAD.

6. Comparar esse modelo como SEIR-QAD compartimental.

Observação: Pode usar uma população menor, como 100 habitantes, para que os cálculos sejam mais rápidos.

Observação: Escolha parâmetros que exibam epidemia.

Referências

• Imperial College COVID-19 Response Team (16/mar/2020), Impact of non-pharmaceutical interventions (NPIs) to reduce COVID- 19 mortality and healthcare demand, 16 de março de 2020. (Free Download)

• Ferguson NM, Cummings DAT, Fraser C, Cajka JC, Cooley PC, Burke DS (2006). Strategies for mitigating an influenza pandemic. Nature 442 (7101), 448–52. (Free Download)

• Ferguson, N. M. et al. (2005), Strategies for containing an emerging influenza pandemic in Southeast Asia, Nature 437, 209-214. (Link para download aqui disponibilizado por CiteSeerX; veja, também, as notas suplementares.)

• Halloran ME, Ferguson NM, Eubank S, et al (2008), Modeling targeted layered containment of an influenza pandemic in the United States. Proc Natl Acad Sci U S A 105 (12), 4639–44 (Free Download). (Veja, também, o suplemento do artigo)

• M. J. Keeling & P. Rohani (2007), Modeling Infectious Diseases in Humans and Animals, Princeton University Press.

• I. Z. Kiss, J. Miller & P. L. Simon (2017), Mathematics of Epidemics on Networks: From Exact to Approximate Models, Springer-Verlag.

Veja a página Imperial College - Covid 19 para mais relatórios, publicações e outras informações do Imperial College.

---

<- Modelos baseados em indivíduos e o relatório do Imperial Colege | Página Inicial | Modelos Individuais Estocásticos Estruturados ->